#ifndef BIO_PH_CELL_INCLUDED
#define BIO_PH_CELL_INCLUDED 1

#include "kernel/solution.hpp"

#include "kernel/ae.hpp"
#include "kernel/nhe.hpp"
#include "kernel/nak.hpp"
#include "kernel/lactate.hpp"

#include "yocto/math/ode/types.hpp"

typedef ode::field<Real>::type DiffEq;

class Cell 
{
public:
	static const size_t NumSol = 4;
	
	static const size_t iS_in  = 1;
	static const size_t iS_out = 2;
	static const size_t idS    = 3;
	static const size_t iS_tmp = 4;
	
	static const size_t NumExtraVar = 1;
	static const size_t iEm         = 1;
	
	
	explicit Cell( lua_State *L );
	virtual ~Cell() throw();
	
	//--------------------------------------------------------------------------
	// num variables = #species(=M) + NumExtra 
	//--------------------------------------------------------------------------
	const size_t M;
	const size_t nvar;
	
	//--------------------------------------------------------------------------
	// Chemical
	//--------------------------------------------------------------------------
	Library          library;
	ChemSys          chemsys;
	Initializer      init_in;
	Initializer      init_out;
	vector<Solution> solutions;
	vector<Real>     variables; //!< nvar, only for Driver
	
	//--------------------------------------------------------------------------
	// Enzymatic
	//--------------------------------------------------------------------------
	EnzymeDB         enzymes;
	AE              *pAE;
	EnzymePtr        EnzAE;
	NHE             *pNHE;
	EnzymePtr        EnzNHE;
	NaK             *pNaK;
	EnzymePtr        EnzNaK;
	Lactate         *pLactate;
	EnzymePtr        EnzLactate;
	
	void             record( const EnzymePtr &enz );
	
	
	//--------------------------------------------------------------------------
	// Geometric
	//--------------------------------------------------------------------------
	const Real      Surface; //!< in microns^2
	const Real      Volume;  //!< in microns^3
	
	//--------------------------------------------------------------------------
	// Electric
	//--------------------------------------------------------------------------
	Real            Em;
	const Real      RT;
	const Real      F_over_RT;
	const Real      Capa; //!< Surface * Cm
	
	//--------------------------------------------------------------------------
	// Formal
	//--------------------------------------------------------------------------
	DiffEq          drvs;
	
	//--------------------------------------------------------------------------
	// Compute net concentrations input
	//--------------------------------------------------------------------------
	void  computeRate( Solution &dS, Real t, Real E, const Solution &S_in );
	
	//--------------------------------------------------------------------------
	// Differential system
	//--------------------------------------------------------------------------
	void  computeDiff( array<Real> &dFdt, Real t, const array<Real> &F );
	mutable Real  dt_ctl; //!< helper for integration
	Real          dt_max; //!< max time step before Newton rescaling
	
	//--------------------------------------------------------------------------
	// For initialisation
	//--------------------------------------------------------------------------
	//Real computeSteadyJ( Real E ); //!< global electric flux at t=0.
	void steadyState();
	
	//--------------------------------------------------------------------------
	// simple h equation
	//--------------------------------------------------------------------------
	void simpleProtonRate( lua_State *L );
	
    Real DeltaQ() const;
    Real DeltaEm() const;
    Real excessConc( const Real Erest) const;
	
private:
	YOCTO_DISABLE_COPY_AND_ASSIGN(Cell);
};

#endif
